Reduced Density Matrix after a Quantum Quench 
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We consider the reduced density matrix (RDM) pA(t) for a finite subsystem A after a global 
quantum quench in the infinite transverse-field fsing chain. It has been recently shown that the 
infinite time limit of pA(t) is described by the RDM pgge.a of a generalized Gibbs ensemble. Here 
we present some details on how to construct this ensemble in terms of local integrals of motion, 
and show its equivalence to the expression in terms of mode occupation numbers widely used in the 
literature. We then address the question, how pA(t) approaches pgge,a as a function of time. To 
that end we introduce a distance on the space of density matrices and show that it approaches zero 
as a universal power-law i~ 3//2 in time. As the RDM completely determines all local observables 
within A, this provides information on the relaxation of correlation functions of local operators. We 
then address the issue, of how well a truncated generalized Gibbs ensemble with a finite number 
of local higher conservation laws describes a given subsystem at late times. We find that taking 
into account only local conservation laws with a range at most comparable to the subsystem size 
provides a good description. However, excluding even a single one of the most local conservation 
laws in general completely spoils this agreement. 



I. INTRODUCTION 



Recent advances in systems of optically trapped ultra- 
cold atomic gases have made it possible to observe the 
nonequilibrium time evolution of isolated many particle 
systems over long time scales^— . A key property of such 
cold atomic clouds is their weak coupling to the environ- 
ment and resulting smallness of external dissipative pro- 
cesses. To a good approximation one is dealing with iso- 
lated quantum mechanical many particle systems, which 
are prepared in generally mixed states, and one is inter- 
ested in the time dependence of observables, in particular 
at late times. The experimental results have stimulated 
theoretical efforts aimed at understanding the principles 
underlying the nonequilibrium dynamics of isolated many 
particle systems. Some of the most basic questions are 
whether observables generally relax to time-independent 
values, and if they do, whether their stationary values 
are described by a statistical ensemble. Other prominent 
issues concern the roles of dimensionality and conserva- 
tions laws. Experiments on trapped 87 Rb atoms^ estab- 
lished that three-dimensional condensates rapidly relax 
to a stationary state characterized by an effective temper- 
ature, whereas constraining the motion of atoms to one 
dimension leads to a much slower relaxation to a non- 
thermal distribution. It was argued that this observed 
difference has its origin in the presence of additional (ap- 
proximate) conservation laws, related to quantum inte- 
grability, in the one dimensional case. Theoretical efforts 
aimed at understanding these and related questions^— 
indicate that in translationally invariant models there are 
at least two basic types of behaviours at late times: sub- 
systems either thermalize^, i.e. are characterized by a 
Gibbs distribution with an effective temperature, or they 
are described by a generalized Gibbs ensemble (GGE) 8 . 
There is evidence that the latter applies to integrable 
models, while the former constitutes the "generic" situ- 
ation. 



A popular protocol for analyzing nonequilibrium evo- 
lution is a so-called quantum quench: here the system 
is originally prepared in the ground state |^o) of some 
local, short ranged Hamiltonian H(ho), where ho is a sys- 
tem parameter such as a magnetic field or an interaction 
strength. At time t = 0, ho is then suddenly "quenched" 
to h, and the subsequent time evolution under the new 
Hamiltonian H(h) is studied. Under this protocol the 
system remains in a pure state \^t) = exp(— iH (h)t)\^o) 
at all times, and as a whole can clearly never be described 
by a Gibbs or generalized Gibbs distribution. This can 
be seen by considering the hermitian operators 



0(«> m ) = \ n )(m\ + \m)(n\, 



(1.1) 



where |n) and \m) are eigenstates of H(h) with energies 
E n and E m respectively. Then the expectation values 

(* t |0(».»O|* t ) = (vK |n)(m|* )e 4 ( £ "-^)* + h.c. (1.2) 

are oscillating in time and never become stationary. A 
useful and intuitive point of view is to focus on local prop- 
erties of a given system in the thermodynamic limit, i.e. 
ask questions only about observables contained in a finite 
subsystem ^l 13 i 31 i 32 . Here the (infinitely large) comple- 
ment A of A can act as an effective bath, and probability 
may freely dissipate from A to A. As a result A may be 
described by a mixed state. Arguably the most precise 
and convenient description of this situation is in terms of 
the reduced density matrix pA{t) of subsystem A. The 
latter is obtained from the density matrix p(t) = \^t)(^t\ 
of the entire system as 



p A (t) = Tr A [p(t)} 



(1.3) 



A central question is then, whether for any finite subsys- 
tem A 



lim p A (t) = /9 B tat,A, 



(1.4) 
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where p s tat,A is a time-independent reduced density ma- 
trix obtained as 



Pstat,A = Tr^ [p 

stat J 



(1.5) 



If (|1.4p holds, then the system evolves towards a station- 
ary state described by the distribution p s ta±- In particu- 
lar (|1.4j) implies that the expectation values of any local 
operator O a acting only within subsystem A is given by 



]im(%\0 A \^ t )=Tr[p stat O a \ 



(1.6) 



An efficient way of investigating whether a given RDM 
approaches a known stationary distribution at late times 
was introduced in Ref. [2(| by considering the operator 
norm || p^(£) — /9 s tat,A Hop- If this approaches zero at late 
times, then the system relaxes locally to the stationary 
distribution p s t a t- Ref- [20| was concerned with the case 
where p s t a t describes a thermal ensemble with a given 
effective temperature, and considered very small subsys- 
tems. Here we arc interested in a quench to a quantum 
integrable model. As we have alluded to before, the sta- 
tionary state for such quenches is believed to be described 
locally by a generalized Gibbs ensemble (for the model 
we consider below this was established in Ref. [32[). More 
precisely, the density matrix of the entire system is ex- 
pected to be of the form 



1 

Pstat = PGGE = 2 C 



(1.7) 



where Z is a normalization 1 , and /„ are local conserved 
quantities, i.e. local operators such that 

[I n ,I m ]=0=[I m ,H(h)]. (1.8) 

The Lagrange mutipliers A„ are fixed by the requirements 



Tr[p. 



(1.9) 



We stress that the GGEs considered in the quench con- 
text are fundamentally different from thermal ensembles, 
because through the specific values of the Lagrange mul- 
tipliers they retain an infinite amount of information 
about the initial state. Above we have stipulated that 
only local (in space) conservation laws I n are to be in- 
cluded in the definition of pgge, but it is in fact a mat- 
ter on ongoing debate whether locality is a necessary or 
even desirable requirement^. In this context a result 
obtained in Ref. 31( is rather illuminating: there it was 
demonstrated for a particular example, the transverse 
field Ising chain, that different statistical ensembles can 
have identical local properties. The two ensembles con- 
sidered were a GGE of the form (|1.7[) . and the so-called 
"pair ensemble" obtained by time averaging the quench 



1 We have in mind regularizing the system by enclosing it in a very 
large but finite volume. 



density matrix p(t). Given that p s ta,t is generally not 
unique, it is clearly desirable to identify the simplest de- 
scription. To that end we introduce truncated generalized 
Gibbs ensembles of the form 



PtGGE 



1 



e "<"o 



(1.10) 



and investigate how well such ensembles describe the sta- 
tionary state for quenches to integrable models. 

The outline of this paper is as follows. In section UU 
we review some relevant results for the transverse field 
Ising chain. Local conservation laws are presented in 
section IIIII and used in sections HVl [Vl I VII to define sev- 
eral classes of generalized Gibbs ensembles. Properties 
of corresponding reduced density matrices are discussed 
in section Ivm In section ["Villi we discuss general proper- 
ties of distances on the space of reduced density matrices 
and introduce the distance used in the remainder of the 
paper. In sections IIX1 IXl IXII and IXIII we present results 
for the distance between quench and generalized Gibbs 
reduced density matrices. We summarize our results in 
section lX!!!! Various technical issues are discussed in four 
appendices. 



II. SOME FACTS ABOUT THE 
TRANSVERSE-FIELD ISING CHAIN (TFIC) 

Here we briefly review some relevant results on the 
TFIC. The latter is an important paradigm for quan- 
tum phase transitions in equilibrium 4 ^ as well as non- 
equilibrium dynamic o 14 ' 24 ' 29 ' 31 ' 32 i 36 ' 49 . In the latter con- 
text experimental realizations range from cold atomic 
gases^ to circuit QED^i. The Hamiltonian of the model 
on a ring is 



H(h) 



3=1 



hat 



(2.1) 



where ct£ +1 = erf. The quantum spins can be mapped to 
(real) Majorana fermions by means of the Jordan- Wigncr 
transformation 




l-i 

au-x = ( ] of, (2.2) 

\3 = 1 



where {a,i,a,j} = 2Sij. In terms of the Majorana 
fermions (|2.2p . the Hamiltonian takes a block-diagonal 
form 



H(h) 



1 + e 



2 

L-l 



1 - e' 



-H 



NS 



H^S/R — ° 2 -? ^ a ' 2 3 + 1 ~ ^ a 2j-l] 

3=1 

—iJa,2L [ha2L—i -F a i] ■ (2.3) 
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Here Af is the number operator 



1. Quenches from the Paramagnetic Phase ho > 1 



(2.4) 



and by construction e l = . crj commutes with 
-Hr.ns- The two blocks Hr and i?NS correspond to 
periodic and antiperiodic boundary conditions on the 
fermions respectively. They can be diagonalized by Bo- 
goliubov transformations 



1 



0*23-1 

a 2 j 



where the Bogoliubov angle 9 p is given by 

h-e* 



Vl + /i 2 — 2hcosp 
The diagonal form of the Hamiltonian is 

#Ns(^) = E Shi?) ( «p"p - o 



(2.5) 



(2.6) 



(2.7) 



pGNS 



where the single-particle energy is given by 



e h {k) = 2J\/l + h 2 -2hcosk. (2.8) 

The ground states of -Hr^nsW are the fermionic vacua 

ap|0;h) RlNS = 0. (2.9) 

Here the momenta are p = where n are even/odd 
integers for R and NS fermions respectively. 



A. Quantum Quenches 

Our quench protocol is as follows: we prepare the sys- 
tem in the ground state \^o) for & n initial value ho of 
the transverse magnetic held. At time t = we instan- 
taneously change the field from ho to h. The state of 
the system at times t > is obtained by evolving with 
respect to the new Hamiltonian H(h) 



(2.10) 



An important quantity is the difference = 9^ — 6^ of 
Bogolibov angles required to diagonalize H (h) and H(ho) 
respectively 



cos Ai 



4J 2 (1 + hh - (h + h ) cos k) 
£h{k)e ha {k) 



(2.11) 



As we are interested in obtaining results in the thermo- 
dynamic limit we have to distinguish between two cases. 



Here the initial state in a large, finite volume is simply 
the NS vacuum 

|*o) = |0;/i o >NS. (2.12) 
The time evolved state can then be written in the form 31 



m = 



aT xp [ z E 



t*nfe)e- 2i ^al p al 



\0;h) 



NS ' 



0<p£NS 

(2.13) 

where |0;/i)ns is the ground state of H^s(h) and M a 
normalization factor. 



2. Quenches from the Ferromagnetic Phase ho < 1 

Here our initial state in a large, finite volume must re- 
flect the spontaneous symmetry breaking of the Z2 spin- 
flip symmetry o~ x ' y — > — a x,y in the thermodynamic limit. 
The appropriate choice is^i 



|0;fe ) R + |0;/i )ns 
V2 



(2.14) 



III. LOCAL CONSERVATION LAWS IN THE 
TFIC 

We consider the one dimensional transverse field Ising 
chain in the thermodynamic limit 



H = -J «+i + h < 



(3.1) 



n— — 00 



Following Ref. [52l | we can construct an infinite number 
of local conservation laws I^r 



= 0, n = 0,l,...,a,/3 = ±, 



(3.2) 



where the Hamiltonian itself is H 
operators 

U n >o 



Uo = 



00 /n—1 \ 

\ e * y. 

^ 00 

-5 E a 3 

j=— 00 



Let us define 



j+n J 



l-i 



U, 



n<0 



\ e < n^k 



j+\n\> 



(3.3) 



j = -oo 



1=1 



and 



Vn>0 



Vn<0 



00 /n—1 \ 

2 e «?n^^' 

j=-OB \l = l / 

1 00 /M-i \ 
"2 E °? II 4«K + |«|- (3-4) 
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In terms of these operators the local conservation laws 
are 



# = -J{U n+l +Ux. n ) + hJ{U n + U. 
J" = J(V n+ i + V-n-i) , n > 0. 



(3.5) 



They are local in the sense that the density of involves 
only spins onn + 2 neighbouring sites. By virtue of their 
locality, the conservation laws can all be expressed in 
terms of Jordan- Wigner fermions 



- Jd2j[a,2j+2n+l + a2j~2n+l], 



-hJa 2 j[a 2 j+2 n—l i&2j — 2n — lj j 



I n -l - -— E a 2j a 2j+2n + a2j-ia2j+2n-l-(3.6) 



We now realize that all conservation laws (|3.6[) are in 
fact quadratic in Majorana fermions! It is then a simple 
matter to diagonalize them simultaneously by means of 
a Bogoliubov transformation (|2.5j) . where on the infinite 
chain the Bogoliubov fermion operators have anticommu- 
tation relations 



{a p ,a{} = 2ir5(p - k). 
The conservation laws take the simple form 



cos(nfc)e(fc)atafc 

27T 



T dp 



2 J sin ((n + l)k)ala k 



(3.7) 



(3.8) 



which furthermore shows that they are even/odd under 
spatial reflections. Interestingly the conservation laws I~ 
do not depend on the transverse field h and are therefore 
shared by the entire one-parameter family of Hamiltoni- 
ans H(h). This seems to be a generic feature of models 
with a free fermion spectrum like the TFIC, see App. [C] 
and Ref. M. 



A. Local conservation laws for periodic boundary 
conditions 

Above we focussed on the bulk contribution to the lo- 
cal conservation laws. For a finite system on a ring there 
are boundary contributions, which can be determined as 
follows. In terms of the Bogoliubov fermions, the conser- 
vation laws for periodic boundary conditions are 



/+ = J2cos(nk)s{k)al 



^ 2 J sin ((n + l)k)a 



k a k 



(3.9) 
(3.10) 



and Fourier transforming back to position space, one ob- 
tains a representation of the conservation laws in terms 
of the Majorana fermions aj for periodic/antipcriodic 
boundary conditions respectively. Finally, inverting the 
Jordan- Wigner transformation gives the desired expres- 
sion in terms of spins. 



IV. GENERALIZED GIBBS ENSEMBLE 

We now define the density matrix of a generalized 
Gibbs ensemble formally by the expression 



PGGE 



\ n=0a=± / 



(4.1) 



where Z is a normalization that ensures TrpcGE = 1. In 
practice we need to regularize (|4.1j) in an asymptotically 
large, finite volume L. 

The Lagrange multipliers are fixed through the re- 
quirements 

h M, (42) 

Using translational invariance we can alternatively work 
with the densities of the conservation laws 



In = E ( J nh 



j'+ra+l j 



to rewrite (|4.2[) as 

Tr [pGGE^n )j,... J+n+l] = (*0 1 ...,j+n+l |*0> 

The solution to this system of equations is 
" dkcos(lk) 



(4.3) 



(4.4) 



A+ = 2-5, 



■arctanh(cos Afc) , 



A," = , 



(4.5) 



where I > and cos is defined in (|2.1ip . In Fig. [T] 
we show for a quench from = 0.1 to h = 0.7. 
The large I behaviour of Eq. (|4.5|) is determined by the 
regions k ~ 0, 7r (where the integrand has a logarithmic 
singularity) and one can show that 



a; 



(-i)' 



(4.6) 



where the sign is + for quenches within the same phase 
and — otherwise. We see that the Lagrange multipliers 
A^ decay rather slowly as a function of n. 



A. GGE in terms of mode occupation numbers 



where the momenta are taken to be either in the R or 
NS sectors. By inverting the Bogoliubov transformation 



In the literature the generalized Gibbs ensemble is of- 
ten constructed from mode occupation numbers n k = 



5 



K 



2.0" 



Here y is an integer and y = 1 (y = oo) corresponds 
to the Gibbs ensemble (GGE). The Lagrange multipliers 
„ are obtained from the requirements 



1.5 



1.0 - 



0.5 



pIgge] 



(*o|(/* 



j,...J+n+l 



(5.2) 

where < n < y. Eqns (|5.2[) are a consequence of 
H ] = and the assumption that the stationary state 
after the quench is described by RDMs based on (|5.ip . 
For transverse field quenches we have X~ y = 0, but the 
other Lagrange multipliers are different from their re- 
spective values in the full GGE, i.e. 



10 



20 



n 



FIG. 1: Parameters X n for a quench within the ordered phase 
from h = 0.1 to h = 0.7. 



a^afc, see e.g4 ' 12 ' 18 ' 53 ' 54 . The latter are non-local (in 
space) as they involve a Fourier transform. We will 
now establish the relation between this and our defini- 
tion (|4.1j) . It follows from (|3.8|) that the density matrix 
can be rewritten in the form 



1 

PGGE = -= exp 



dk t 
2^ lk k0lk 



(4.7) 



where 



7fe = ^2 K £ h{k) cos(fcn) 

n=0 



2J\: 



(k(n + l)). (4.8) 



This establishes the fact that the GGE density matrix 
can be constructed either from the local conservation 
laws (|3.6[) . or from the mode occupation numbers n^. 
This relationship generalizes to interacting integrable 
models, where the appropriate GGE can be formulated 
either in terms of the local integrals of motion generated 
from the transfer matrix, or from the mode occupation 
numbers n a {k) = Z\{k)Z a {k) : where Z a (k) are Faddeev- 
Zamolodchikov operator o 55 ' 56 . 



\+ / \ + 



(5.3) 



We note that the correlation matrix of can be com- 
puted efficiently using FFT algorithms. This is in con- 
trast to the case of theories with non-trivial scattering 
matrices, for which it is extremely difficult to reconstruct 
the Lagrange multipliers from the conservation laws^i. 



VI. DEFECTIVE GENERALIZED GIBBS 
ENSEMBLES 

It is instructive to consider a second type of truncated 
GGE, where we retain an infinite, but incomplete set of 
integrals of motion. Such "defective" GGEs will allow us 
to ascertain the role of a particular local conservation law. 
We define the truncated defective GGE as the density 
matrix (q < y) 



,(+q),V 

"tdGGE 



1 



■ exp 



T+1 

n,(+q),y n 1 



(6.1) 



in which the Lagrange multipliers / + > arc fixed by 
the system (|4.2[) with n < y, n ^ q, and we have used 
that the Lagrange multipliers A~ ^ +q ^ y must vanish as 
a consequence of reflection symmetry around the origin. 
We then define the defective GGE as the limit y — > oo of 
truncated defective GGEs: 



PdGGE 



(6.2) 



V. TRUNCATED GENERALIZED GIBBS 
ENSEMBLES 

In order to assess the importance of the various con- 
served quantities, it is useful to define ensembles that in- 
terpolate between the Gibbs distribution and the GGE. 
We define particular such truncated GGEs as follows. 
Given that the densities of the conservation laws 1^ in- 
volve n + 2 neighbouring sites, it is natural to retain only 
the "most local" conservation laws, i.e. 



V n=0 a=± 



(5.1) 



In order to solve the system of equations (|4.2p for the 
defective GGE it is useful to work in the mode occupation 
number representation (|4.7[) . which reads 



/'dGGE 



1 

exp 

Z(+q) 



dk 



(6.3) 



where the Lagrange multipliers 7^ 
set of equations 



(+q) 



are subject to the 



dfc 
2^ 

dk 
2~H 



tanh 



tanh 



cos Aj 



(+?) 



cos A 



e(k) cos(nfc) = , 
sin((n + l)fc) = 0.(6.4) 
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Guided by the fact that cos(nfc) and sin((n + 1)A;) form 
an orthonormal set of functions on [— 7r, 7r], we look for a 
solution of the form 



tanh 



7^ 



= cos Afe 



+ cos(gfc) 
iq e{k) 



(6.5) 



where rej" is a yet to be determined constant. We note 
that the value of k+ affects the expectation values of local 
operators. For some cases k+ can be easily determined 
as follows. Given that | tanh(cc)| < 1, (|6.5[) implies that 



cos A i. — K 



, cos(qfc) 



< 1 



Vfc. 



q s(k) 
Setting k = 0, n then gives 

\2J(h - l)sgn(ft - 1) - K+| < 2J\h - 1 
\2J(h + l)-(-l) g K+\ < 2J(h+l) 



(6.6) 



(6.7) 



Eqs (|6.7p allow us to identify cases, in which k+ = 0: 

1. odd q and quenches within the same phase; 

2. even q and quenches across the critical point. 

Importantly, k+ = implies that P^q E = Pgge, i-e. the 
defective GGE is identical to the full GGE. This "GGE 
reconstruction" is a peculiarity of free-fcrmion models 
and can be traced back to the existence of conserva- 
tion laws independent of the quench parameter, see also 
Ref. [E3. 

For general quenches and values of q, is determined 
by Eq. (|6.2p . We find that it takes the value correspond- 
ing to the maximal entanglement entropy (as shown in 
Fig- although the entanglement entropy may be non- 
stationary under a variation of the excluded integral of 
motion (see Appendix |D] for further details). 



Jordan- Wigncr transformation (|2.2[) . The nonlocality of 
the transformation (|2.2p has important consequences for 
RDMs. First and foremost, if the spins are not adjacent, 
the map from spin to fermionic degrees of freedom does 
not have a simple reduction to the subspace of the Hilbert 
space formed by sites {x\,--- ,xi}, because of Jordan- 
Wigncr strings stretching between site a 59 i 60 . However, 
the RDM of a block of adjacent spins can be mapped 
one-to-one on a block of adjacent fermion a 61 i 62 , provided 
that the first site of the block coincides with site 1, i.e. 
the origin of Jordan- Wigner strings. Then (|7.1|) can be 
represented in the form 



P£ 



'21 



(7.2) 



with Hj = 0, 1. An important quantity in what follows is 
the correlation matrix T 



Tij = Tr [pajai] 



1 < i,j < 21 



(7.3) 



In the cases of interest to us, the correlation matrix is of 
block- Toeplitz form 



r = 



r 



r_i 
r„ 



(7.4) 



where 



n dk 



-ilk 



-f(k) g(k) 
-g(-k) f(k) 



(7.5) 



The TFIC Hamiltonian exhibits a Z2 symmetry 



VII. REDUCED DENSITY MATRICES IN THE 
TRANSVERSE FIELD ISING CHAIN 

In this section we summarize some basic features of 
RDMs in the TFIC. We note that most of the follow- 
ing discussion generalizes straightforwardly to other spin 
chains with free fermionic spectra such as the quantum 
XY model. Our starting point is a density matrix p de- 
scribing the entire system, which we take to be of size 
L with periodic boundary conditions. We are interested 
in the limit L — > 00, but it is convenient to start with a 
large, finite chain. The RDM of a subsystem consisting 
of I spins 1/2 at sites Xi, i = 1, . . . ,£ can be expressed in 
the form 



P{xi,-" ,x e } 



-y 



Tr[p< •■•<]<••■<, (7.1 



where a, = 0, x, y, z and er° = I. The quantum 
spins are mapped to (real) Majorana fermions by the 



aj — > -aj . (7.6) 

If the density p is invariant under the transformation 
(|7.6[) we have Tr [p aj] = 0, and as Wick's theorem applies 
to the Jordan- Wigner fermions we can express (|7.2j) as a 
Gaussian^ 



Pi 



1 

Z 



-e 4 



I £ m „ a m W mn a n 



(7.7) 



where Z ensures that Trp£ = 1 and W is a skewsymmet- 
ric, 2€-by-2i! Hermitian matrix related to T by 



tanh — 
2 



(7.8) 



We now turn to the three particular cases of interest, 
namely those where p in (|7.1[) is a thermal density ma- 
trix, a GGE density matrix, or the density matrix after 
a global quantum quench of the transverse field in the 
TFIC. 
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A. Thermal Density Matrix 



D. Defective GGE Density Matrix 



On a very large ring, the Hamiltonian has a block diag- 
onal structure, see section|TTJ The thermal density matrix 
is a function of the Hamiltonian and therefore inherits the 
same block structure 



Pp 



1 



1 



(7.9) 



2 Z-r, 2 Zns 

It follows from this that only even operators have non- 
vanishing expectation values, i.e. 

Tr (ppO) ^ -> [e**^, O] = 0. (7.10) 

In the thermodynamic limit the difference between ex- 
pectation values of local operators with respect to the 
R and NS sectors tends to zero, so that we may work 
exclusively in e.g. the R sector. The resulting RDM of 
a contiguous block of spins is then Gaussian (|7.7[) , (|7.8[) 
with 



(I>L=Tr 



Z 



R 



— Si 



It can be written in the form (17.41) with 



m = o 



k tanh 



(7.11) 



(7.12) 



B. GGE Density Matrix 



It was shown in Ref. [32j (see also [24j|31|) that the 
RDM of the generalized Gibbs ensemble (|4.1 j) .(|4.7 |) is 
Gaussian and can be expressed in the form (|7.7j) . (|7.i 
The correlation matrix is given by 



CTggeL - ^Tr 



j ajdi 



(7.13) 



n 3 z 

It can be written in the form (|7.4[) with 

/(fc) = , g(k) = -ie Wk tanhf— 

Here the 7fc's are related to the A^'s by (|4.8p and the 
Bogoliubov angle 9k is given in (|2.6[) . 



(7.14) 



C. Truncated GGE Density Matrix 

The correlation matrix of the truncated GGE defined 
in Sec. |V]is given by 



(r& E ).. = ^Tr 



y-i 

ex p(EE[ A Vn])«i 



- Si. 



(7.15) 



n=0 a=± 

It can be written in the form (|7.4[) with 

/(fc) - , g{k) = -ie Wfc tanh(P v -i(a>s(i0)e(*)) . 

(7.16) 

Here P y —i(x) is a polynomial of order y — 1, which is 
computed numerically. 



In Sec. ED wc defined the defective GGE P^qqe &s 
le ensemble that lacks in the 
correlation matrix is given by 



the ensemble that lacks in the conservation law J+. Its 



(+<?) 
dGGE 



r(d) 



Tr 



exp 



T1=Q 



(7.17) 

It can be written in the form (|7.4[) with /(fc) = and (cf. 
Eq. GH) 



tanh 



Ik 



+ cos(qk) 



(7.18) 



where k+ is computed numerically maximizing the en- 
tanglement entropy, which selects A^~, +g s = whenever 
it is allowed. We note that the Fourier transform of 
Eq. (|7.18p . which is required to compute the correlation 
matrix (|7.5[) , can be easily expressed in terms of the GGE 
correlators; for \£\ < q we have 



T dk 



-ilk 



g(k) = 



n dk 



-ilk 



3GGE(fc)- 



^sgn(log/i)^- 1 e-l lo s /l l < ? 



(7.19) 



Since is a bounded function of q (cf. Eq. (|6.7[0 . at 
fixed £ the fermionic correlators approach the GGE ones 
at least exponentially fast in q. 



E. Quench Density Matrix 

At zero temperature the ground state phase diagram 
of the TFIC exhibits ferromagnetic (h < 1) and param- 
agnetic (h > 1) phases, separated by a quantum critical 
point. In the ferromagnetic phase the Z2 symmetry of 
the Hamiltonian is broken spontaneously. As we will see, 
this symmetry breaking has important effects on the time 
evolution of the density matrix. 



1. Quenches originating in the paramagnetic phase 

Here at t > the full quench density matrix is 

p{t) = \m t ) (* t | , (7.20) 



where the state \^>t) is given in (|2.13|) . As a result of 
the squeezed-state form of \*f?t), Wick's theorem applies 
to averages calculated with respect to p(t), and RDMs 
are Gaussians of the form (|7.7[) . (|7.8[) . with correlation 
matrix 



r(i) - ns (0|e lffNst a jaj e- iffNst |0) rTO - Si 



'ns "n 



(7.21) 



8 



This is of the form (17.41) with 



by the Frobcnius norm 2 



cos Afe — i sin Afc cos(2eh{k)t) 



f(k) = sinAfcsin(2e? l (fc)i), 
where e l6k is given by (|2.6[) . 



(7.22) 



2. Quenches originating in the ferromagnetic phase 



Given the initial state ()2.14|) . the post-quench density 
matrix of the full system is 



p(t) = |* t )(* t | , 

e- iH ^\0;h ) K + e 



0; ho) 



NS 



V2 



Importantly, RDMs are no longer Gaussian in this case. 
We will discuss how to cope with this complication in 
section IXII It is knowr*^ that in the stationary state 
RDMs are Gaussian with a correlation matrix equal to 
the t -> oo limit of ([7T2T]) . 



VIII. DISTANCES ON THE SPACE OF RDMS 

In the following we focus on RDMs for finite subsys- 
tems of lattice models with a finite dimensional Hilbcrt 
space at each site. In this case the RDMs are finite di- 
mensional matrices and a simple way to define a distance 
between two density matrices is by means of a matrix 
norm 



da{p,p') =|| p-p' \\a 



(8.1) 



Here the index a labels different matrix norms. As we 
are dealing with finite matrices, all norms are equivalent 
in the sense that 



Cab || P \\a<\\ P \\b< C h a II P IU 



(8.2) 



where c a b and q, a arc positive numbers that depend on 
the matrix dimension but are independent of p. One 
consequence of (|8.2p is that if the distance between two 
matrices approaches zero when some external parameter 
p is tuned to a value p, the dependence of the distance 
onp-p is almost independent of the norm. On the other 
hand, the dependence on matrix dimension is in general 
very different for different norms. This is important for 
our purposes, because the matrix dimension is related 
to the size of the subsystem under consideration, and it 
is principally desirable to be able to compare distances 
between different sizes. 

From a technical point of view, the distance induced 



U\\i 



Tr [A^ A] 



(8.3) 



is generally the easiest to calculate. On the other hand, it 
has the drawback that the physical interpretation of the 
distance is less transparent than for some other norms. 
For instance, given two density matrices p and p' , a very 
natural question is how different expectation values of lo- 
cal observables are in the two ensembles. We now discuss 
this question for the particular case of spin- 1/2 quantum 
spin chains. Here the most important local observables 
are products of Pauli matrices. These are particular cases 
of involutions O 2 = I, for which the following inequality 
holds 



\Tr[(p-p')6]\ 



< o- 



P- P 



(7.23) Here 



Tr [VZ4t] 



(8.4) 



.5) 



is the trace norm. From Eq. (|8.4p it is evident that the 
trace distance provides an upper bound for the difference 
between the expectation values of observables in the two 
states: if | p — p' ||i< e, then the expectation values 
of all (local) observables will agree in the two ensembles 
within accuracy e. In terms of the Frobcnius distance we 
have instead (here we use that the local Hilbert space is 
two-dimensional) 

|Tr[(p-p')0]|<||p-p'||i<2 f / 2 ||p-p'|| F . (8.6) 

On the other hand, we have 

|Tr [{p - p')6] | < |Tr [ P 6] \ + |Tr [p'6] \ < 2, (8.7) 

where in the last step we have used that for involutions 

6 

|Tr [pO] | < £ | {Ph I = trV^P = 1- (8-8) 



Combining (|8.7[) and (|8.6j) we see that as long as || p — 
p' ll-F~ 2 1 ~ e / 2 , the Frobenius distance does not provide 
useful information about expectation values. It is shown 
in Appendix [5] that for sufficiently large £ this is always 
the case. 

A second problem with using the Frobenius norm as a 
distance is that the norms of RDMs at late times after 
a quantum quench, as well as the norms of RDMs de- 
scribing Gibbs or generalized Gibbs ensembles, generally 
are exponentially small in the subsystem size. Given the 
upper bound derived in Appendix [5] 



P 



P\\% 



P' 



2 



(8.9) 



2 It is also known as the "Hilbert-Schmidt norm", but we prefer 
to call it "Frobenius norm" to emphasize that we are considering 
finite subsystems. 
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this implies that in these cases of interest | p — p' \\p is 
exponentially small in subsystem size. This shows that 
the Frobenius norm itself is not a convenient measure for 
the distance between two density matrices. The same 
problem occurs for the operator norm ||A|| op = \Amaxi 
where A max is the largest eigenvalue of A. This norm 
was used for example in Ref^ to analyze the relaxation 
properties of small subsystems after a quench into a non- 
intcgrable model. Indeed we have 



II P~ P' \\op<\\ P Hop + || P 



I op 



(8.10) 



and the maximal eigenvalues of RDMs for large subsys- 
tems are generally exponentially small in subsystem size. 



A. Definition of the Distance 

In order to circumvent the problem described above, 
we define our "distance" 3 on the space of RDMs as 



Wp-p'Wf 



VTp 



P'\\ F 



.11) 



1. An upper bound: Using the upper bound derived in 
(|A1[) we see that 



V{ P ,p') < 1. 



.12) 



2. A lower bound: A lower bound for T>(p, p') can be es- 
tablished by means of the triangle inequality | \p— p' \ \p > 
\\p\\f — \W\\f I- Using that the Frobenius norm of a 
RDM is related to the second Renyi entropy by 



5 2 = -logTr[p 2 ] =-log||p||| 
we find that 



3.13) 



B. The Distance between two Thermal Ensembles 



In order to establish a benchmark for (|8.11[) , it is useful 
to consider the distance between the RDMs of two ther- 
mal ensembles at slightly different inverse temperatures 
P and j3' (but the same Hamiltonian). Then 



V{p!},pp>) 



1 



1 1 dp ft 1 1 

WppWf V2 



.16) 



For a sufficiently large subsystem (and a local Hamil- 
tonian), the first factor on the right hand side can be 
expressed as 



ll gPg 112 

II dp I If 

INI! 



\\ Pp ((h) -h)\\ 2 f 



\pp\\f 



Tr 



Trpj 
(((H), -Hf) 



20' 



(8.17) 



where (0)p = Tr^O). For a large subsystem, this is 
proportional to the square of its size, and hence 



hgrllf; 



t 



(8.18) 



We conclude that the distance between two thermal 
RDMs on a subsystem of size £, and (3 sa /?' is 



V(pp, P(s .)<xl\p-p'\ 



(8.19) 



As expected this is proportional to the difference in tem- 
peratures, but there is also a factor of I, The latter is 
important if one is interested in comparing the distance 
between two ensembles for different subsystem sizes. 



\P-P'\\> 



S 2 

cxp [ - — 



S' 2 

ex P [ -— 



(8.14) 



The Distance between two GGEs 



This provides the desired lower bound 

| e -S 2 /2_ e -S^/2| 



?.15) 



We note that the bound (|8.15[) is independent of subsys- 
tem size I as long as the second Renyi entropies of the 
two ensembles differ at least by a constant (when viewed 
as functions of £). 



3 We have not proven that T>(p,p') obeys the triangle inequality 
as this is not essential for our purposes. 



The above discussion carries over to the case of two 
generalized Gibbs ensembles (|4.1[) , with slightly different 
values of Lagrange multipliers A^. The leading contri- 
bution to the distance is given by 



V (pgge,Pgge) 



II dPGGE II . 

y-y ll a r m H F l 



|PGGE||F 



V2 



(8.20) 



A calculation similar to the thermal case shows that for 
large subsystem size 



I dpc 



3A° 



IPggeIIf 



•21) 
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D. Information on observables contained in the 
distance 

Let us consider the situation where the distance be- 
tween two reduced density matrices p\ and pi defined 
on an interval of length I becomes small, and denote the 
corresponding averages of local operators on said interval 
by 



(0) a = Tr [ Pa O] , a = 1,2. 



5.22) 



By expanding the density matrices in a complete basis of 
Hermitian involutions we can show that 



v(0 0) - EoM^mi 



5.23) 



Defining an average 



P{0) 



(0)l + {0) 2 2 



E Q (Q) 2 i + (Q)V 



.24) 



we can express the distance (|8.11[) as 



V( Pl ,p 2 ) = [(R(0)) 



1/2 



Here 



R{0) 



(0)\ + (0)1 



.25) 



.26) 



is the relative difference between the ensembles described 
by pi and p 2 , respectively. This implies that T>( P \,p2) 
measures the mean relative difference of the expectation 
values of all local operators, averaged with respect to the 
probability distribution (|8.24|) . 



E. The Distance between two Gaussian Density 
Matrices 



The distance (|8.1ip between two Gaussian RDMs p[T] 
and p[V] can be expressed in terms of their correlation 
matrices following Ref. (59[. Given the definition of the 
distance 



T>(p,p') = 



Tr(p 2 + p' 2 - 2pp>) 



5.27) 



/Tr(p 2 ) +Ti-(p' 2 ) 
we require tractable expressions for the quantities 

Tr(p[r>[r']). (8.28) 



This is achieved in two steps. First, we note that the 
product of two Gaussian RDMs (|7.7p is itself Gaussian^ 



exp 



^7 X! 11 j cxp fi 



CLidj 



'■J 



exp 



En^' 



WW 



y\ij Ojidj 



(8.29) 



'•j 



This can be seen by expanding the left hand side of (|8.29[) 
by means of the Baker-Campbcll-Hausdorff formula in 
terms of multiple commutators, and then observing that 
the commutator of quadratic operators is quadratic and 
gives rise to a commutator between the matrices W that 
characterize the 2-forms 

[\ E W ij a i a h \ E WijOiOj] = i E| 

i,j i,j i,j 

(8.30) 

Second, using results for the second Rcnyi entropy^, one 
can relate the Frobcnius norm of a Gaussian RDM to the 
correlation matrix by 



:[[w, 



.31) 



Tr[p[r] 2 ] = (det|l±^^ 
Combining (|8.31j) and (|8.29[) one can then show^ that 

i + rf 



{r,r} = Tr[p[r] P [r]] = det 



32) 



Here we have used that Tr [P1P2] > 0, which is a con- 
sequence of density matrices being positive semidefinite 
operators. Finally, substituting (|8.32[) into (|8.27p . we 
obtain the following result for the distance between two 
Gaussian RDMs 



W],/>[f]) 



2{r,f} 



{r,r} + {r,r} 



.33) 



Given that the correlation matrices are only 2£ dimen- 
sional (with £ the subsystem size) , (|8.33[) provides a very 
efficient way of computing distances for large subsystem 



sizes. 



IX. SINGLE-SITE SUBSYSTEM 

It is instructive to consider the time evolution of the 
RDM describing a single-site subsystem in some detail. 
In this case the RDM of site 1 can be expressed in the 
form 



pi(t) = ^+m(t) -o x 



(9.1) 



where fh(t) is the magnetization per site at time t after 
the quench, i.e. 



m«(t) = i(* t |<7? |* t ) 



(9.2) 



11 



The RDM of the generalized Gibbs ensemble describing 
the stationary state is 



Pgge.i = - + m z Btat al 



where 



in' 



71 dk 



(9.3) 



(9.4) 



Finally, the RDM of the thermal ensemble described by 
pp, whose inverse temperature f3 is fixed by the require- 
ment 

lim ~(*o|#0)|*o> = lim yTr [ Pl3 H(h)}, (9.5) 



is given by 



.1 = 2+ m /3 a i 



(9.6) 



Here the transverse magnetization per site is 

m% = ^ge^tanh(^). (9.7) 

A. Quenches originating in the paramagnetic phase 

Here the Z 2 symmetry enforces 

m x (t) = m y (t) = 0. (9.8) 

The z-component of the magnetization per site is 

/*■ (fa 
— e l9k [cosA fc -isinA fe cos(2e fc i)] . (9.9) 

For late times we may evaluate the integral by means of 
a stationary phase approximation, which gives 



m z (t) ~ m z stat + 



c(t) 



where 



c(t) = 



(Jt) 3 / 2 ' 

(h - h ) cos(4Ji|l -h\- tt/4) 
8\h -l\y/Tr\h-l\ 

(h - h ) cos(4Ji|l + h\ + tt/4) 
8\h + lWn\h + l\ 



(9.10) 



(9.11) 



The distance between p\{t) and the generalized Gibbs 
RDM at late times then decays to zero like a power-law 
with exponent 3/2 



^(/Ol(i),PGGE,l) = 



y/2\m*(t) - m z 



Vl + 2(m^)) 2 +2(m^ GE ) 2 



l + 4(m* tot )< 



On the other hand, the distance between pi(t) and the 
thermal RDM approaches a constant at late times 



y/2\m z {t) - m 



) Jl + 2(m z (t)) 2 +2(m^) : 



^l + 2(mf tat ) 2 +2(m^' 



+ 0((Jt)~t). (9.13) 



B. Quenches originating in the ferromagnetic phase 

Here all three components of the magnetization per 
site are non-zero. The component along the transverse 
field direction is again given by (|9.9|) . while the late-time 
asymptotics of m x (t) has been calculated in^i 



1 



™ x (t) = -VCf^e^o 



(9.14) 



Here Cp F is a known amplitude and e' k = dE jl • Finally, 
the Hcisenberg equation of motion for af (t) relates the y 
and x components 



m y (t) = 



1 dm x {t) 
2Jh dt ' 



(9.15) 



Importantly, m x ' v (t) exhibit exponential decay in time. 
In contrast, m z (t) again decays like a power law with 
exponent 3/2 and therefore will dominate the late time 
behaviour. Hence at sufficiently late times, the distances 
of pi(t) to GGE and thermal RDMs arc again given by 
(|9.13[) and f|9.12[) respectively. So for a single site subsys- 
tem the spontaneous symmetry breaking only modifies 
the intermediate time behaviour of the distances. As we 
will see, this holds true also for larger subsystems. 



X. LARGER SUBSYSTEMS FOR QUENCHES 
FROM THE PARAMAGNETIC PHASE 

For quenches with ho > 1 and in the thermodynamic 
limit, we determine the distance between the quench 
RDM and that of an appropriate thermal or generalized 
Gibbs ensemble by means of relation (|8.33j) . The correla- 
tion matrices for all cases are of the form (|7.4[) , (|7.5|l with 
elements given in (|7.12p . (|7. 14[) and (|7.22[) respectively. 
For a subsystem of size I this requires the calculation 
of determinants of 21 x 21 matrices, which is done nu- 
merically. Results for a quench from ho = 1.2 to h = 3 
and subsystem sizes I = 10, 20, 30 ... , 150 are shown in 
Figs [2] and [3] We see that the distance between quench 
and Gibbs RDMs tends to a ^-dependent constant at late 
times. This establishes that subsystems do not thermal- 
ize. On the other hand, as can be seen from Fig. [31 
at sufficiently late times the distance between pi{t) and 
Pgge,£ decays to zero in a universal power-law fashion 



T^(pe(t),PGGE,i) 



Jt»i 



k(£){Jt)~ 3 / 2 



(10.1) 



12 



A O V «. 



10 
150 




V 



(GGE) 




FIG. 2: Normalized distance £> Gibbs = V(j>t(t),p$) after a 
quench within the paramagnetic phase for subsystem sizes I = 
10, 20, . . . , 150. As £ increases, the color fades from brown to 
green, the symbols become smaller and the curves narrower. 
At late times the distances tend to constants depending on 
subsystem size. 



The quality of the fit (fTUTj) is shown in Fig. gj The 




10 100 1000 

t h =1.2-»h=3 



FIG. 4: Distance £> GGE = V{pt(t), pf ) after a quench 
within the paramagnetic phase for two representative values 
I — 70, 150. We used the same notations of Fig. [3] The black 
dashed curves are best fits to the form T> = at~ 3 ^ 2 . 



6000 - 



4000 - 

t y 

2000 - ^ 

if 

- v <3a£> - 0.01 



100 200 300 

£ h =1.2-+h=3 

FIG. 5: Dependence of time on subsystem size at fixed dis- 
tance V(pi(t), pf GE ) = 0.01 for the same parameters as in 
Fig. [3] The dashed curve is the best fit to the functional form 
t = a + b£ 4/3 . 



FIG. 3: Normalized distance £> GGE = V(p e (t), pf GE ) af- 
ter a quench within the paramagnetic phase for subsystem 
sizes £ = 10, 20, . . . , 150. As £ increases, the color fades from 
brown to green, the symbols become smaller and the curves 
narrower. At late times D(pe(t), pf GE ) tends to zero in a 
universal power-law fashion oc (Jt)~ 3 ^ 2 . 



large-^ asymptotics of the function k(£) can be inferred 
as follows. On surfaces with constant, small T> the time 
scales as t ~ £ 4 ^ 3 as is shown in Fig. [5J This in turn 
implies that 



k(£) 



(10.2) 



A. Relaxation time 

We may extract a relaxation time from the behaviour 
of the distance, by using the connection to averaged dif- 
ferences in the expectation values of local operators es- 
tablished in subsection I VIII rJ] The distance can be writ- 
ten as 



V(p, 



GGE, 



,p e (t)) = ([R(0)i 



1/2 



where 



R{0)= l<°>*-<°>GGB 



y/(0) 2 t + (0) 



2 

GGE 



(10.3) 



(10.4) 
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and the bar denotes the average (|8.24[) . Using that 

Rip) < \j[R{O f = V( PGGEJ , p t (t)), (10.5) 

and then substituting the asymptotic behaviour (jlO.ip . 
(|10.2[) into the right hand side, we obtain 

R(0) < £ 2 t~ 3/2 . (10.6) 

Bounding the right hand side by a (small) constant, we 
obtain a time scale t* ms associated with the relaxation of 
the average relative error with respect to the distribution 

C^ 4/3 - (10-7) 

It is not simple to identify the observables that give signif- 
icant contribution to the average, since it depends both 
on their "multiplicity" in the subsystem (produced by 
translational invariance and other symmetries) and on 
the expectation values. We note that the relaxation time 
t* ms is very different from the time scales identified in 
Ref. [32| in the time evolution of the two point functions 
of spin operators for quenches within the paramagnetic 
phase. 




1 10 100 

£=10 / h =1.2-»h=3 



FIG. 6: Distance V {v) = T>{pt{t), p$ GB e ) at fixed length 
£ = 10 between quench and truncated GGE reduced density 
matrices for y = 1, 2, 4, 8, 16 and a quench within the param- 
agnetic phase. Here y is the maximal range of the densities 
of local conservation laws included in the definition of the 
ensemble. As the number of conservation laws is increased, 
the time window, in which the distance decays as t~ i ^ 2 , in- 
creases. At very late times all distances with finite y saturate 
to nonzero values. 



B. Distance from Truncated Generalized Gibbs 
Ensembles 

Having established that the distance between quench 
and GGE reduced density matrices tends to zero as a 
universal power law at late times, a natural question is, 
how close the quench RDM is to the truncated GGEs 
(|5.ip . which retain only finite numbers of conservation 
laws. A representative example for a quench within the 
paramagnetic phase is shown in Fig. [5] We see that at 
sufficiently late times, the distances converge to constant 
values. However, increasing the range (and number) of 
conservation laws, the values of these plateaux decrease, 
signalling that retaining more conservation laws gives 
better descriptions. In an intermediate time window, the 
extent of which grows with y, the distance decays in a 
universal t~ 3 l 2 power-law fashion. In Fig. [7] we consider 
the distance 

(10.8) 

between the RDMs of the truncated and full generalized 
Gibbs ensembles as a function of the parameter y. For a 
given subsystem size £, this corresponds to plotting the 
values of the plateaux seen in Fig. [6] against the corre- 
sponding values of y. The distance is seen to start de- 
caying exponentially as a function of y as soon as y > £ . 
There are two main conclusions of the above analysis: 

1. Including more local conservation laws improves 
the description of the stationary state. 




O v v o 

J 1 1 1 1 I 1 1 1 1 Lli 1 A LOj 1 1 1 I 1 l_j i_l 

10 20 30 40 5C 



y h =1.2-»h=3 

FIG. 7: Distance T>^ — D{pGGE,t, Puzge e) between the 
GGE and the truncated GGEs obtained by imposing local 
conservation laws with densities involving at most y + 1 con- 
secutive sites. The quench is from ho = 1.2 to h = 3 and 
the subsystem size ranges from I = 5 to I — 50. Colors 
and sizes change gradually as a function of the size I. For 
y > £, the distance starts decaying exponentially in y, with 
an ^-independent decay constant. 



2. The description of the stationary state improves 
rapidly, once the range y + 1 of all conservation 
laws not included in the truncated GGE exceeds 
the subsystem size £. 
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C. Distance from defective Generalized Gibbs 
Ensembles 

We now turn to the role played by particular local con- 
servation laws. We find that the distance between quench 
and defective GGE reduced density matrices for a given 
quench and subsystem size tends to a constant at late 
times, i.e. 

^^ftW^L,^*' 1 ' (10-9) 

The dependence of this asymptotic value on the subsys- 
tem size £ and the integer q is shown in Fig.[8]for a quench 
within the paramagnetic phase. We see that V^ +q) ex- 



vy 



1x10 
i(+q) 



□ 5 

50 



40 60 

h n =1.2->h=3 



FIG. 8: Distance Ut^ +q ^ = "D{pgge,£, P^jge l) f° r a quench 
within the paramagnetic phase, for subsystem lengths I — 
5, 10, ... , 50. The excluded conservation law is 7+ with even q. 
Colors and sizes change gradually as a function of the length. 
When y > £, the distance starts decaying exponentially with 
a decay length given by Eq. (|10.10p . 



hibits an exponential decay in q as soon as q > I. This is 
similar to the behaviour observed in the truncated GGE 
case. The decay length can be calculated from the large-*? 
asymptotics of Eq. (|7.19|) . By series expanding Eq. (|8.33[) 



to second order in T 



(+9) 
dGGE 



Tgge we obtain 



Numerically we find that k+ ~ l/<7 2 - 



(10.10) 



1. "GGE Reconstruction" 

In section IVII we discussed the issue that, for certain 
quenches and omitted conservation laws /+, the corre- 
sponding defective GGE is identical to the full general- 
ized Gibbs ensemble. We now return to this point. In 
Fig. [5] we consider the truncated, defective GGE for a 
quench across the critical point from ho = 2 to h = 0.5 
for a subsystem of length £ = 5. We plot the distance 



between the reduced density matrices of the appropriate 
GGE and the truncated, defective GGE with y integrals 
of motion, where J+ (q < y) has been excluded, i.e. 



T)d(+q),y T)( „ J.+l),V \ 



(10.11) 



As discussed in section IVII for even q we expect this 
distance to approach zero, when the number y of con- 
servation laws goes to infinity. This behaviour is clearly 
observed in Fig. |H1 This implies that the corresponding 
conservation laws do not affect averages of local opera- 
tors. As discussed before, this is a particular feature of 
free theories, where H(ho) and H{h) generically share 
certain local conservation laws. 
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FIG. 9: Distance V^ +qhv = V( Pg ge,i , P^GGE.t) for a 
quench across the critical point between the GGE and the 
defective truncated GGE RDMs for a subsystem of 5 consec- 
utive spins, as a function of the number y of retained conser- 
vation laws. Each symbol corresponds to a different excluded 
conservation law Iq (the legend indicates the value of q). The 
distance approaches zero for even q, whereas it remains finite 
for odd q, in agreement with the discussion of section fVTl The 
lines are the distances from the corresponding defective gen- 
eralized Gibbs ensemble pJjgge with maximal entanglement 
entropy. 



2. More local conservation laws are more important 

On the other hand, for odd q we find that Vii +q) ' y 
approaches constant values when y becomes large. This 
value agrees with the distance between the GGE and the 
defective GGE with maximal entanglement entropy (we 
stress that for the considered quench the defective GGE 
does not always correspond to a stationary point of the 
entanglement entropy under a variation of the excluded 
integral of motion, as shown in Fig. [20] of Appendix [D]) . 

The fact that V < ^ 3 +q ^ ,v tends to a constant at large y 
shows that retaining an infinite number of local conser- 
vation laws while excluding one of them is insufficient for 
describing the stationary state. By comparing distances 
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for different values of q we observe that for a given value 
of y, T>'^ +q ' > ' v decreases as a function of q. This implies 
the more local the conservation law, the more important 
it is for describing the stationary state. 



XI. QUENCHES FROM THE 
FERROMAGNETIC PHASE: EFFECTS OF 
SPONTANEOUS SYMMETRY BREAKING. 

We now turn to quenches originating in the ferromag- 
netic phase, i.e. ho < 1. In this case, the time evolved 
initial state is given by 



l^t>R+ l^t>NS 



iH a t 



V2 

« t \0;h ) a a = R,NS, (11.1) 



where |0; /io)r/ns are ^ e g roun d states of the Hamilto- 
nian H(h ) with pcriodic/antiperiodic boundary condi- 
tions. In order to analyze reduced density matrices after 
a quantum quench from the ferromagnetic phase we will 
make use of the following facts. 

a) The fcrmion parity e™ N = flj G ] <EH> is full y 
factorized in space. 



b) The states \tpt)n and IV^ns are eigenstates of e 
with eigenvalues 1 and —1 respectively. 



c) The difference between the expectation values of local 
operators in the states IV^r and I^Ons tends to zero 
in the thermodynamic limit. 

d) The RDMs Tr A [|?/> t ) aa (V>t|], a = R,NS, where A is a 
single interval and A its complement, arc Gaussian. 

Property [b| allows us to express the full density matrix 
in the form (c/. Eq. ([TTjO 



P(t) 



;2Re[ NS (^iaiVt> R ]a}, 



O, 



(11.2) 



where Z ensures that Ti[p(t)] = 1 and {O e } U {0 } is a 
complete set of Hermitian involutions with the property 



0. 



{e i7Tjv ,O a } = 0. 



(11.3) 



We will refer to O e / as even and odd operators respec- 
tively. The main difference between even and odd opera- 
tors, is that the latter are not local in terms of fermions: 
a Jordan- Wigncr string is attached to them. We are in- 
terested in the RDM of a block A of £ contiguous spins, 
which is obtained by tracing out the degrees of freedom 
outside A 



(11.4) 



A convenient representation for pi is obtained by restrict- 
ing the sums in Eq. (|11.2[) to involutions that act as the 
identity operator outside the interval A, i.e. 



O -> {A) ® I^ 



(11.5) 



where the superscript (A) indicates that the operators 
act on the Hilbert space over all sites in A. As a result of 
property [3$, fermion parity has a simple restriction onto 
the interval A 



leA 



i ■ 



(11.6) 



and can be used to subdivide operators into even 

and odd ones 

[e i ^ A ,O{ A ^]=0, {e^,OW}=0 . (11.7) 

This then implies that we can decompose the RDMs of 
(|11.2[) into even and odd parts as well 



Pi = PLe + Pl,o 



(11.8) 



In the thermodynamic limit we then may employ prop- 
erty [g) to obtain the following expressions 

PiAt) = ^Y. Rc [^t\0 \4' t ) R ]0 . (11.9) 



Importantly, the even part pe t e(t) is Gaussian (|7.7j) by 
virtue of property |d|, and has the same structure as 
RDMs for quenches originating in the paramagnetic 
phase. On the other hand, the odd part pi i0 has 
its origin in the spontaneous breaking of the Z 2 sym- 
metry. The commutation relations pi.7[) imply that 
Tr [pg }0 (t)pf a ] = for any Gaussian density matrix pf 3 -, 
because the latter is by construction even. As a result 
the odd part pi, of the RDM enters the distance from a 
Gaussian state only through its norm 



V( Pe (t),p^) 




Pt,e(t)-Pt* \\ 2 F + \\PtAt) 



12 
IF 



\\f + \\ PLo \\ F + \\ P^Wf 

(11.10) 

We will be interested in the cases where pf & describe 
Gibbs or (truncated) generalized Gibbs ensembles. The 
Frobenius norms \\pt^ e {t) — pf a \\F, \\pe,e\\F and Hp^Hf 
can be efficiently evaluated by means of Eq. (|8.33|) . What 
remains in order to determine the distance (|11.10[) is a 
method for calcuating the Frobenius norm || pg_ \\p. 
This is a somewhat involved technical problem, which is 
addressed in Sec. IXI Al and Appendix [Bj The basic idea 
is to utilize a cluster decomposition theorem at any finite 
time after the quench, see also Ref. (36j . 
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FIG. 10: Geometry of the composite system Au n used in 
calculating || pi, ||f, where pi is the RDM of subsystem A. 
The single site at position n is separated from A by a block 
B of length r. 



A. || pt t0 \\f from cluster decomposition 

The main difficulty in calculating the Frobenius norm 
of pi^ is that the latter is not Gaussian. The idea is 
therefore to obtain p£ jD as a reduction of a Gaussian op- 
erator. To that end, we consider a composite system 
C = A U 11 consisting of our subsystem A and a single 
site at position n, which is separated from A by a block 
B of length r, see Fig. [ITJ 

The even part of the RDM pc(t) can be expanded in 
a complete basis of Hermitian involutions O e / Q as 



pc At) = ^[^Mo^ 



EE 

O a—x,y 



{0 a a n )0 a a n , (11.11) 



where (• • ■) = (* t | ■ • ■ |\P f ) « R (ip t \ ■ ■ ■ IV'^R, since both 
O e (j^ and (j" in (jll.lip are even operators with re- 
spect to fermion parity. In the limit of large separation 
r, we may use the cluster decomposition principle to sim- 
plify the expectation values 

(O e <) {0e ) «> , 

(0„O ^ ( 0o )«). (11.12) 

This then leads to the following relation between RDMs 
in the limit or large separation 

lim p C , e (t) = Pe,e( t )®PlAt) + Pt,o( t )®P-LAt), (1H3) 
r— >oo 

where pi is the RDM of site n. The piece of interest to 
us is 



pe,o <X> pi,o 



lim 



1 



O a a—x,y 



(0 aZ)0 a%. (11.14) 



In the next step we move from spins to Major ana 
fermions by means of the Jordan- Wigner transformation 

(El 



Pt,o®f>l,o= lim TjTXT Y. KAie l ^ B )A e™ M *a 



ittAJ b „a 
n i 



(11.15) 

where A are odd products of Majorana fermions acting 
on sites within A. Importantly, the fcrmionic expres- 
sion f| 1 1 . 1 5[) depends on the configuration of Majoranas 



in subsystem B through the Jordan- Wigner string oper- 
ator. The right hand side of (| 1 1 . 1 5[) can be cast in the 
form 



Pi,o ® pl,c 



lim ( e ^B) e ^B? — In^n ^ (n^g) 



where p is a normalized, Gaussian operator (|7.7[) acting 
on the Hilbert space over sites A U n 



P 



(11.17) 



In writing (|11.15[) we are assuming (e I7rA/ ^) ^ 0. The fact 
that p is Gaussian is a consequence of the particular form 
of \ipt) R (which is the analog of (|2.13j) in the R sector) 
and Mb being quadratic in fermions. The odd part of 
the single-site RDM is of the form 



Pl At)=™ X (t)<+™ y (t)° V n 



(11.18) 



and hence 

12 



[piAt)? = (KW1 2 + K(t)] 2 )l 2 = ml(t)I 2 . (11.19) 

Here the late-time behaviour of m x ' v (t) are give n by 
(Pn%|) and (|9~T5j) respectively, and following Ref. [3l[ they 
can be easily calculated numerically for all times. Com- 
bining (|11.19[) and (|11.16|) we obtain 



vkNe 



Pi,o || F 



= lim 



>l 



r^oo y/2\m±(t)\ 



9 - KWn 



lim — ; r— - 

r->oo 2|mx(t)| 



Tr[t*-(v> 9 y 



(11.20) 



Since both p and er*p<7* are Gaussian, their moments can 
be written in terms of their respective correlation matri- 
ces 



Gij = Tr [fa j a,] - 5y , 
Oij = Tr [a^a^a 3 ai] - 5 i: 



(11.21) 



We note that the correlation matrices are related by Q = 
PnQPni with P n the diagonal matrix that changes the 
sign of the last 2-by-2 block (Id is the d x d identity) 

Pn = he®(-h)- (11-22) 

Using (|8.32p we have 

Tr[p 2 ] ={Q,G}, Tr[«p) 2 ] ={Q,Q}. (11.23) 

A slight complication arises because p is not positive 
semidefinite. To account for this we must use the more 
general definition of {T, T'} as the product of the eigen- 
values of (1 + rr')/2 with halved degeneracy^. We may 
then recast pi.20p in the form 



II PLo \\f= lim 



2\m±(t)\ 



{G,G}-{g,g}. (n.24) 
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r h =0.2-h=0.8 

FIG. 11: Difference || 8p e , \\f = \\ pi,o \\f — lim || pi :0 \\ F as 

r— >oo 

a function of the separation r for a quench from ho — 0.2 to 
h — 0.8. We see that for r > 2n max i the difference becomes 
exponentially small in r/£, where £ is the correlation length 
in the initial state. 



II Pe.o \\f is significantly more costly, and in practice in- 
volves determinants of at most 2(£ + 2v max t + £5) x 2(£ + 
2v max t + £S) matrices, as discussed above. 

Results for a quench from /io = l/3to/i = 2/3 and sub- 
system sizes £ = 10, 20, 30 ... , 150 are shown in Figs Q21 
and [131 We see that the distance between quench and 
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While formally correct, (|11.24[) is not suitable for numer- 
ical computations, because at large distances (e Z7rJ ^ B ) be- 
comes very close to zero. A more convenient expression 
derived in Appendix [B] is 



lim 



det(l 2 ^ © 2r © h + iT A uBun) 



2 1 + e / 2 \m±(t)\ 

(11.25) 

Here Taubuu is the correlation matrix of the interval 
A U B U n and is given by (|773|) . ([73]) . ([T22]) . In order 
to utilize pi. 251) we in principle have to consider infinite 
separations r and hence infinitely large matrices. 

Crucially, in practice a finite separation r > 2u max t, 
where v max = maxfc£/j(fc) is the maximal propagation ve- 
locity, is sufficient to recover the r — > oo limit up to cor- 
rections that are exponentially small in r/£. Here £ is the 
correlation length in the initial state. A representative 
example is shown in Fig. 1111 In practice, using a finite 
r > 2v max t + £6 with 6 ~ 20 provides an efficient way 
for calculating ||/)f, (t)||i? and then by means of (|11.10[) 
distances T>(pe(t), pf a ) for quenches originating in the 
ferromagnetic phase. 



FIG. 12: Distance £> Gibbs = V(pi(t),pf) after a quench 
within the ferromagnetic phase for subsystem sizes £ = 
10, 20, . . . , 150. As £ increases, the color fades from brown to 
green, the symbols become smaller and the curves narrower. 
At late times the distances tend to constants depending on 
subsystem size. 



Gibbs RDMs tends to a ^-dependent constant at late 
times. On the other hand, as shown in Fig. [T31 at suffi- 
ciently late times the distance between pe(t) and pgge,£ 
decays to zero in a universal power-law fashion 



V{pi{t),PGGV,l) 



k(£)(jty 3/2 + 



(11.26) 



The large- £ asymptotics of the function k(£) can be in- 
ferred in the same way as for quenches within the param- 
agnetic phase. On surfaces with constant, small T>, time 
scales as t ~ £ 4 ^ 3 as is shown in Fig. Q3J which implies 
that 



k(£) 



(11.27) 



B. Results for quenches from the ferromagnetic 
phase 

For quenches with ho < 1 and in the thermodynamic 
limit, we determine the distance between the quench 
RDM and that of an appropriate thermal or general- 
ized Gibbs ensemble by means of relations ([ll.lOp and 
(|11.25p . The correlation matrices for all cases are of the 
form (|7T4|) . ((731) with elements given in ([7X2]) . (|7TT4|) and 
(|7.22j) respectively. For a subsystem of size £ most terms 
require the calculation of determinants of 2£ x 2£ matri- 
ces, which is easily done numerically. The evaluation of 



We conclude that the late time behaviour of the dis- 
tance between quench and generalized Gibbs RDMs is 
the same as for quenches within the paramagnetic phase. 
The mean relaxation time i* ms is therefore again given 
by p0.7p . Interestingly this coincides with the result ob- 
tained in Ref. (32j for the relaxation of the order param- 
eter two-point function after quenches within the ferro- 
magnetic phase. The effects of the spontaneous symme- 
try breaking are important only at short and intermedi- 
ate times. It is shown in the inset of Fig. [13] that there is 
a time window, in which the odd part of the RDM gives 
the dominant contribution to the distance, which decays 
exponentially. 
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FIG. 13: Distance E> GGE = T>(p e {t), pf GE ), after a quench 
within the ferromagnetic phase for the subsystem lengths I = 
10, 20, . . . , 150. We used the same notations of Fig. [3] The 
behavior is almost the same as that shown in Figs [3] and [5] 
but the effect of the spontaneous magnetization is visible at 
intermediate times, when the distance decays exponentially 
(inset). 




odd part pn Q of the density matrix. The relative impor- 
tance of pf jD for large I can be estimated by considering 
the von Neumann entropy of subsystem A 



S v n[pi] = Tr[/9 f \n(p e )} 



(11.28) 



We recall that the von Neumann entropy after a global 
quench grows linearly in time until the Fermi time tp = 
£/(2v max ), and then saturates to a value proportional to 
the subsystem size l 64 i 65 . Using the commutation rela- 
tions (|11.8[) we see that the even part pi_ e can be ex- 
pressed in terms of the full RDM pi as follows 



Pi 



e nN ' A p l e i ^ A 



(11.29) 

Since for any set of density matrices pi the von Neumann 
entropy satisfies^ (A.^ > 0, J^. \ = 1) 

A* log Xi < S vN XiPi] - ^ [pi] < , 

i i i 

(11.30) 

the following bounds on the von Neumann entropy of 
subsystem A hold 



SvN[pe,e] - log 2 < S vN [p e ] < SvN^e]- 



(11.31) 



This demonstrates that at any time after the quench the 
symmetry breaking contribution to the von Neumann en- 
tropy will be at most log 2. Given that for large subsys- 
tems the von Neumann entropy at late times is propor- 
tional to I, we conclude that the relative contribution of 
the odd part of the RDM will be important only for small 
subsystem sizes. 



FIG. 14: The time vs. the subsystem length at fixed distance 
V(p e {t),pf GE ) = 0.01 (black solid line of the left plot). The 
dashed curve is t = a + bl 4 ^ 3 , with a and b obtained by fitting 
the numerical data. The filled region shows the effect of the 
spontaneous magnetization. 



C. Magnitude of the contribution due to p( j0 

The effects of the spontaneously broken Z2 symmetry 
in the initial state make themselves felt through the Z2- 



1. A conjecture for \\ p£ t0 \\f in the limit of large £ and Jt 
We now consider the space-time scaling limits 

£,Jt^ 00, — fixed. (11.32) 

We observe that in this limit our numerical results for 
quenches within the ferromagnetic phase are in excellent 
agreement with the following relation 



f n dk 

log || p t , {t) log II p t , e {t) \\ F + / — log (cos Afc) max {0,2s' k t - £ + 0(£°, t )} . (11.33) 

I 

Here we have highlighted the asymptotic nature of the portant corrections will arise. Since log | |p^ e (i)| \p is pro- 
relation and indicated by O(£ ,t°), where the most im- 
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portional to the Renyi entropy S2 (c/. Eq. (|8.13I0 . wc 
may use the known results^ 4 - on the asymptotics of the 
latter 

log || peAt) I|f= -S2/2 M 

* dfc bg l + cos 2 A fc m . n(24 ^ t) + Q{£0> t0) _ 



2tt 



(11.34) 



Combining (|11.34p and (|11.33|) provides a conjecture for 
the asymptotic behaviour of || pi_ Q \p. This conjecture is 
compared to numerical results in Fig. [15] The agreement 
is clearly quite good. 





FIG. 15: The ratio R = !^-°^ll F after a quench within the 
ferromagnetic phase for subsystem lengths I = 10, 20. The 
lines correspond to the analytic expression (|11.33l) . where we 
have included a correction 0{l ) by shifting i — > I — 1.2. The 
inset presents the same data on a logarithmic scale. 



FIG. 16: Distance E> GGE = V(p e (t), pf GE ), after a quench 
from ferromagnetic phase to the paramagnetic phase for the 
subsystem lengths I = 10, 20, . . . , 150. The conventions are 
the same as in Fig. [3] 
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XII. QUENCHES ACROSS THE CRITICAL 
POINT 



FIG. 17: Distance £> GGE = T>(p e (t), pf GE ), after a quench 
from ferromagnetic phase to the paramagnetic phase for the 
small subsystems I = 1, 2, 3, 4. 



We now turn to quenches across the critical point. 
These are of particular interes t 14 ' 31 i 42 . In Fig. [16] wc 
plot the distance between quench and GGE reduced den- 
sity matrices for a quench from the ferromagnetic phase 
(ho = 1/2) to the paramagnetic phase (h = 3/2). The 
15 data sets displayed correspond to subsystem sizes be- 
tween £ = 10 and i = 150. We find that the distance 
£>gge _ D(pt(t),pf GE ) again decays in a universal i~ 3/2 
power law. In Fig. |TB|wc consider the same quench, but 
focus on very small subsystem sizes ^ = 1,2, 3, 4. We ob- 
serve that the distance displays an oscillatory behaviour 
on top of a power-law decay in time. This is in agreement 
with the analytic results discussed in section HXB I for the 
I = 1 case. Increasing the subsystem size leads to a rapid 
suppression of the amplitude of the oscillations. 

In Figs [TH] and [TH] wc consider the reverse quenches, 
i.e. starting at ho = 3/2 in the paramagnetic phase, and 



quenching to h = 1/2 in the ferromagnetic phase. The 
behaviour of the distances is very similar to what we 
found for the quench from h$ = 1/2 to h = 3/2: at late 
times the distance decays as a t~ 3 / 2 power law, and for 
small subsystem sizes we observe oscillatory behaviour 
on top of this decay. 



XIII. SUMMARY AND CONCLUSIONS 

In this work we have considered the evolution of re- 
duced density matrices after a quantum quench in the 
transverse field Ising chain. The main result of our work 
is to demonstrate that 

lim p e (t) = p G GE,e, (13.1) 
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t h =3/2-»h=1/2 

FIG. 18: Distance T>(pi(t), pf GE ), after a quench from para- 
magnetic phase to the ferromagnetic phase for the subsystem 
lengths t = 10, 20, . . . , 150. The conventions are the same as 
in Fig. |3] 




10 20 30 40 50 



t h =3/2->h=1/2 

FIG. 19: Distance T>(pt(t), pf GE ), after a quench from para- 
magnetic phase to the ferromagnetic phase for the small sub- 
systems I = 1, 2, 3, 4. 



where pe{t) is the reduced density matrix of a subsystem 
consisting of I adjacent spins after a quench of the trans- 
verse field, and pgge,£ is the reduced density matrix of an 
appropriately defined generalized Gibbs ensemble. The 
derivation of (|13.1j) is based on defining an appropriate 
distance T>(p,p') on the space of reduced density ma- 
trices, and then establishing that the distance between 
quench and GGE reduced density matrices approaches 
zero at late times. For our particular choice of distance 
we found that at late times this distance approaches zero 
as a universal power law in time 

2?(Mi),PGGE,£)~£~ 3/2 ■ (13.2) 

We have presented a detailed construction of PGGE,e 
in terms of the local (in space) integrals of motion 1^ of 
the TFIC. The densities of these conservation laws in- 
volve only spins on n + 2 consecutive sites. We proved 



that these local conservation laws are related in a lin- 
ear fashion to the occupation numbers of the Bogoliubov 
fcrmions that diagonalize the Hamiltonian of the TFIC. 
This linear relation establishes the equivalence of our con- 
struction of the GGE to the one frequenctly used in the 
literature, which is based on mode occupation numbers. 

We then have addressed the question, which of the 
conservation laws are most important for obtaining an 
accurate description of the stationary limit lim^oo pe(t) 
of the quench RDM. To that end we introduced (defec- 
tive) truncated generalized Gibbs ensembes, which are 
missing some of the local conservation laws. We found 
that the more local the conservation laws (i.e. the fewer 
consecutive spins their densities involve), the more im- 
portant they are for describing the stationary state of 
a given subsystem. Loosely speaking we observed that 
in order to obtain a good description of the stationary 
state RDM of a subsystem of size £, we need to retain all 
local conservation laws, whose densities involve at most 
pa t + n neighbouring spins, where n is a constant de- 
pending on ho and h. Leaving out "highly local" conser- 
vation laws generally provides a very poor description of 
the stationary state. To the best of our knowledge this 
is the first such demonstration of a connection between 
locality of conservation laws and their importance in the 
GGE context. 

Our work raises a number of issues. First and foremost 
is the dependence of the results obtained on the precise 
definition of the distance on the space of reduced density 
matrices. We have argued, that the "best" distance is 
the one based on the trace norm, because it provides the 
most direct and precise information on the time evolution 
of local observables. Unfortunately this distance is much 
harder to handle analytically. It would however be very 
interesting to implement it in purely numerical studies 
using iTEBD or related algorithms. 
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Appendix A: Inequalities involving the Frobenius 
norm of RDMs for spin-1/2 Quantum Spin Chains 

In this appendix wc provide lower and upper bounds 
for the Frobenius norm of the difference of two reduced 
density matrices || p— p' \\f in a translationally invariant 
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system. An upper bound is obtained as follows 

\\p-p'\\% = T,[p 2 +p' 2 -2pp'] 

= \\p\\ F + \\p'\\ F -2Tr(pp>) 

< Up 111 + Up' III- (Al) 

Here we have used that both p and p' are positive 
semidefinite and hence 



— A m i n Trp' — A m j n > 0, 



(A2) 



where < A m ; n < \j are the eigenvalues of p. To derive 
a lower bound we start by expressing the RDM of a block 
of I spins in a spin- 1 chain in the form 

\ i Y J ^\p°T---°?\°T---o7 , (A3) 



H ~ 2 l 

{a,} 

where a, = 0, x, y, z with er° = I, and p is the density 
matrix of the full system; pi is only function of the length 
because of translational invariance. By singling out the 
term with ai = 0, we can express this in the form 



Pi 



3 

-1 <£> I X at a e 



(A4) 



ai—l 



where pt-\ is the RDM of the block consisting of sites 
!,...,£-!. We also write the RDM of the I th spin 



I 3 

pi = , + E Tr [ s pZi\ °\ 



and observe that 

|| Pl -p' 1 f F =2j^{Tr[nZ 1 ]y 



(A5) 



(A6) 



at — 1 



Here we have defined fi^ = 8p^ - 6p'"J v Using (|A4 
we have 



- p'i IIf 



P*-l - P<>_1 \\F 



2E ll^li 



> 



P£-i - Pe-i \\ F 



3 ^[O^]) 2 



E 



2 f-2 



Pi-i-ti-i IIf , II Pi - Pi IIf 



2 *-i 



(A7) 



where we have used that for N x N matrices M we have 
NTtM 2 > (TtM) 2 in the second step, and (|A"6l) in the 
last. Iterating Eq. (|A7[) £ — 1 times we obtain 

(A8) 



\\Pi-p' t \\ F > 2'-H\\p 1 -p' 1 \\ F . 

This implies that for sufficiently large subsystem size £, 
the distance || p£ — p' t \\ F will generally be larger than 



Appendix B: Derivation of Eq. (|11.25p 

Our starting point is Eq. (|11.24[) . i.e. 



Pl,o II F= hm 



r^oo 2|mj.(i)| 
Our task is to evaluate 

{G,G} and (e 



{G,g}-{g,g}. (bi) 



{G,G}, (B2) 



where G = PnGPn and P n is the diagonal involution 
defined in (|11.22[) . We recall that {Q,Qi} denotes the 
product of the eigenvalues of (I + QQ\)/2 with halved 
degeneracy (the eigenvalues of QQi are always double 
degenerate^ for antisymmetric matrices Q and Qi). The 
correlation matrix G defined in Eq. (|11.21[) turns out to 
be^ the Schur complement of the block matrix of the 
matrix r^uBun, 

G = TAUn — Fauu, B Fb, AUn ■ (B3) 
t B 

Here i r 2 denotes the matrix, whose rows and columns 
are associated with spatial regions R\ and R2 respec- 
tively, e.g. 



Tij+2t 0<i<2£ 

Fi+2r,j+2£ i > 21 . 



(B4) 



Here £ is the size of subsystem A, while r is the distance 
between A and site n, see Fig. [TU1 

The first step is to find determinant representations for 
{Q, Q} and {Q, PQP}, where P is a generic symmetric 
involution (P 2 = I and P t = P). 

We first consider {Q, Q}. Since Q is antisymmetric, its 
eigenvalues come in pairs ±q 

= dct \Q - ql\ = det |Q* - ql\ = det | - Q - ql\ . (B5) 

Both eigenvalues ±q give rise to the same eigenvalue 1+q 2 
of I + Q 2 , and hence 



{Q,Q} = \ 



1 + q 2 



(B6) 



<?>o 



Here the product is over all positive eigenvalues of Q. 
Using that 

det |I+<Q| = H(l + iq) H(l -iq) = l[(l + q 2 ) , (B7) 

<?>0 q > q>0 

it follows that 

{Q,Q} = det \(l + iQ)/V2\. (B8) 
Next we consider {Q, PQP}. The matrix 



p-2 



1 e™/ 4 ! + c~ i7f / 4 P 



V2 



(B9) 



22 



satisfies (P 2 ) 2 = P and (P 2 )' = P 2 - Since we have 



I + QPQP = (P 5 )"^ 1 + (P*QP>) 2 )P* 



(BIO) 



the eigenvalues of l + QPQP and 1+{P 2 QP 2 ) 2 coincide. 
Therefore 



{Q,PQP} = {P 2 QP 2 ,P 2 QP 2 } = det 



I + iPzQP? 



V2 



where in the last step we used Eq. (|B8|) . Since (P 2 ) 2 
and P 2 = I, (|B11| can be rewritten in the form 



(Bll) 



P 



{Q, PQP} = det \P\ det \(P + iQ)/V2\ . (B12) 



Using (|B8|) and (|B12[) , we can reexpress the quantities in 
(lB2l) as follows: 



(e 4 ^) 2 {Q, G} = det |*T B | dot |(I + iS)/V2\ 

(e^f{g,g} 



det \iT B \ det \(P n +i 



G)/V2\ 



(B13) 



Here we have use that the expectation value of the string 
operator in region B is related to the correlation ma- 
trix T B by (e t7rAAf3 ) 2 = det|irs|. A remaining problem 
is that lim r _ ) . o0 det \iT B \ = 0, which precludes a numeri- 
cal evaluation of (|B1[) on the basis of expressions (|B13D . 
This complication is overcome as follows. We recall the 
expression of the determinant of a block matrix 



det 



Mu M12 
M 2 i M 22 



det|M 22 |dct 



Mu - M 12 M^M 21 



(B14) 

We then substitute (|B3[) into (|B13p . and identify 

{<?, <?} as the de- 



2 i+i ^Mb} 1 fagy and 2^+1 
tcrminants of the matrices 

i + ir^un ii\4 UniB \ / p„ + ir^un il^u^s' 
«r B .Aun «r B / \ «r SiJ 4 Un iT B 

(B15) 

respectively. Rearranging some of the rows and columns 
we obtain 



{G,G} 



det \he © 2r © I 2 + ir^usun | 
2^+i 



/ l7 rAA B \2 rr , det \ht © 2r © (~I 2 ) + iT AuBUn \ 



2 e+1 

(B16) 

The representations (|B16[) are suitable for numerical cal- 
culations even in the limit of large r. There is one further 
simplification: in the limit r — > oo we have 

lim (e l ^f{g,g} = - lim (e MA ^) 2 {£, £}. (B17) 

r— >oo i — >oc 



To see this, we expand the determinants in (|B16I) with 
respect to the last 2-by-2 block (from here on we omit 
the subscript in r^usiM, i.e. T = T A uBun) 

det \l 2l © 2r © I 2 + iT\ + det |I 2f © 2r © (-I 2 ) + iT\ = 

2 det|r + zl 2 ^ © 2r+2 | - 2dct|r Au s + ihe ©0 2r | . 

(B18) 



Using properties of the correlation matrix one could show 
that the determinants on the second line approach zero 
in the limit of large distance. For the sake of simplic- 
ity we propose a different proof, which is based on the 
assumption that the limit 



lim dot \Taub + ihe © 2 J 



(B19) 



exists: we demonstrate that the limit cannot be infi- 
nite, so the expression in Eq. (|B18|) docs tend to zero as 
r -> oo. To this end we consider the {21 + 2r) x {21 + 2r) 
correlation matrix G of a generic Gaussian density ma- 
trix, and show that the determinant det \G + il 2 e © 2r | 
has an upper bound independent of r. Hence it cannot 
diverge in the limit r — > oo. Our proof is based on the 
following facts: 

(a.) || G \\op< 1, and hence || G 2 \\ op < 1 and || G + iI 2 z® 
2l - ||o P <||G|| op +l<2; 

(b.) G + il 2 t © 2r cannot have more than 21 eigenvalues 
with absolute values exceeding 1. 

Property ju) is a consequence of G being the correlation 
matrix of a positive semidefinite Gaussian. Property (|b]) 
can be proved as follows: Let w a normalized vector with 
Wi = for any i < 21. Then 

vft(G + il 2 i © 2r ) f (G + Hu © 2r )w = w^G 2 w < 1 , 

(B20) 

where the inequality follows from property (Jaj) . If there 
were more than 21 eigenvalues A of G + il 2 e © 2r with 
modulus larger than 1, wc could find a linear combination 
W = J^i c iVi of the corresponding normalized eigenvec- 
tors Vi with the property Wi = for any i < 2£; this 
leads to a contradiction with (|B20[) since 

c *4 (G + H 2e © 2r ) t (G + il 2t © 2r ) c :>Vj 

i j 
i i 

This completes the proof of property (JbJ ) . 
Properties (jaj) and JbJ) imply that 



det|G + z/ 2 £©0 2r || < 2 



2/ 



(B22) 



which establishes that the determinants in (|B18|) remain 
finite in the limit r — > oo. Concomitantly the expression 
in Eq. (|B18[) approaches zero as r — > oo. This establishes 
(|B17[) . Putting everything together we see that (|B16p 
can be written as 



pe,o \\f 



= lim 



det l 2 i © 2r © I 2 + iT AuBUn 



2 J 2-+ 1 \m ± (t)\ 

(B23) 

which is Eq. (|11.25j) . 

We stress that our assumption reagrding the 
limit (|B19[) is equivalent to the existence of the limit 
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in ()B23[) . From a numerical point of view, this can be 
inferred from the scaling analysis of 



det he © 2r © h + i^AuBim 



2I+ 1 ! 



-(*)! 



(B24) 



which is still required to check the cluster decomposition 
hypothesis (see Fig. [TTj). 

The magnetization |mj_(i)| can be computed writing a 
self-consistent equation for Eq. (|B23I) in the case I = 1: 
From Eq. (|11.19[) we have 



|| Pho ||f.= >/2|ro±(f)|, 
which together with Eq. (|B23[) gives 



4mi(i) = lim Wdet(I 2 © 2r © I 2 + ir JuBL 



(B25) 



(B26) 



Appendix C: Conservation laws in spin models with 
free fermion spectra 

In this appendix we present a simple construction of 
the bulk contribution to local conservation laws of the 
TFIC on the infinite line. Our method readily general- 
izes to other models with free fermionic spectrum such as 
the XY chain. Ignoring boundary conditions, we can use 
the Jordan- Wigner transformation to express the Hamil- 
tonian as a quadratic form in Majorana fermions 



H = — ^ a{Hi n a n 



(CI) 



£,71 

Here H is a skewsymmetric block-circulant matrix 

yo yi ■■■ y L -i 
y L -x y I 



H = 



Jo 



(C2) 



where y n = — jj_„ are 2-by-2 matrices. In Fourier space 
we have 



I^ e ^(y fe) ,.,, 



(C3) 



k=l 



where {Y^)j n = —(Y_k) n j. One can show that a complete 
set of local conservation laws is obtained by taking 



1 x - 

l,n 



(C4) 



From Eq. (|8.30j) we see that [H,I r ] — if and only if 
[H,I r ] = 0. Similarly one has [I r ,I r i] — if and only if 
[I r ,I r r] = 0. Hence the problem of constructing conser- 
vation laws is equivalent to determining an appropriate 



set of mutually commuting matrices that commute with 
H. Because the projectors on the eigenvectors of blocks 
circulant matrices are block circulant matrices, we seek 
X r in block-circulant form 



yir) yir) 



y\ 



(r) 



x, 



(r) 



(C5) 



Imposing [H,Xr] = and \Z r ,X r i\ = we obtain the 
conditions 



M v(r'h 



Vfc. 



(C6) 



where Y^ r ' is the Fourier transform (|C3[) of y^ . In the 
quantum Ising model arc 2-by-2 traceless matrices, so 
Eq. (|C6[) has the simple solution 



Y, 



(r) 



(C7) 



where to 



(r) 



-ijj_ u and 



'I 



Fourier transforming 



back to position space we have 



y. 



1 L 



2i"k n (r) T 



fc=i 



1 ^ ^ 2™k ( r ) 

fe=i 



(C8) 



We define the 'range' of a local conservation as the maxi- 
mal number of neighbouring spins involved in its density 
minus one. By construction, the range is equal to the 
maximal |n| such that y^ is nonzero (cf. Eqs (|C1|) . 
([C2])). For the TFIC one finds that Y n = for |n| > 1, 
and concomitantly the range of the Hamiltonian is r# = 
1. It is straightforward to identify the conservation laws 
with ranges < r + 1: from Eq. (|C8|) they are such that 



r+l 

Wfc = c~ sin(nfc) , 

n=l 



r+l— m 



(Ik 



! = C « C0S ( nfc ) 



n=0 



(C9) 

They can be divided in two classes: one with q k = 0, 
which we denote by J - , and one with u>k = 0, which we 
denote by / + . Finally, a complete set of conservation 
laws is given by 



i + : y 



ij2e 2 -^ n co S (rk)Y k 

fe=i 



y-,(r) 



2J V— ^ 2lYik 



(CIO) 



L 



sin((r + l)fc)I 



fe=i 



These are exactly the conservation laws reported in 
Eq. (l2~T3l . 

We note that the conservation laws 7~ are independent 
of the system details, and can be found in any noninter- 
actin g m odel with a block circulant structure (see also 
Ref. [58|). Indeed they are originated from the trivial 
solution of Eq. (|C6|) , namely the identity. 
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Appendix D: Peculiar aspects of defective GGEs 

In this appendix we discuss some properties of the de- 
fective generalized Gibbs ensembles defined in ScctionfVT] 
We start by recalling the standard variational approach 
for deriving statistical ensembles in quantum mechanics. 
One generally seeks the density matrix that maximizes 
the entropy under a given set of constraints on indepen- 
dent, additive conservation laws I 3 

STr [-p log p - Xp - £ Xjljp] = . (Dl) 

3 

The solution of (|D1[) is of the form p cx expQ~V ^jlj): 
which shows that the ensemble is a function only of the 
conservation laws appearing in Eq. (|D1I) . 

We now consider the density matrix after a quench. 
All the ensembles defined in the main text are compati- 
ble with the principle of maximal entanglement entropy, 
and the GGE, the truncated GGE, and the truncated 
defective GGE can be obtained (a posteriori) by means 
of the variational approach (|D1I) . 

Some complications arise when we consider defective 
GGEs, in which we exclude a single integral of motion. 
From Eq. (|6.5p we find that the entanglement entropy 
density a^ GE< " +9 ^ of the defective GGE /0^qq E 4 is given 
by 

dGGE(+q) r Ak Tj( A + COs((jfc)\ 

a vN »=y o _Jj(cosA fc -<-^), (D2) 

where H(x) = — LLe log — i^logi^. By writing 
the defective GGE as in Eq. (|4.2[) . one can easily show 

^dGGEf + q) 

that avN + is the Lagrange multiplier associated to 

Ofcq 

the conservation law J+ (c/. Eq. (|4.5[) L if the maximum 
of the entanglement entropy is not at the boundaries of 

the domain of fc+ then the equation — s ^t T = has 

" Ofcq 

a solution, and jO^GGE can obtained from Eq. (|D1I) . 
In the absence of peculiar constraints, one would expect 



the maximum to be generally a stationary point of the 
entanglement entropy. However, quenches in translation- 
ally invariant noninteracting models are very special since 
the initial state is a simultaneous eigenstate of an infinite 
number of local conservation laws. This substantially re- 
duces the degrees of freedom, and can result in an ex- 
ceptionally small domain for fc+ (which may not include 
a stationary point). In Fig. 1201 we show this paradoxical 
behaviour for the same set of parameters used in Fig. [9] 
Besides the pathological cases of even q, in which the 
curves collapse to the point k+ = 0, the effect of the 
reduction of degrees of freedom is reflected in the "trun- 
cated" shape of the curves for q ^ 1, which turn out to 




q 2 K + q h =2-»h=1/2 



FIG. 20: The difference of entanglement entropy densities 
Aa V N = <y^% GE ^ +q ^ — 0"{fjv E as a function of the parameter 
Kq for the same quench shown in Fig. [5] (the legend indicates 
the value of q). The points have the maximal entropy and 
correspond to the lines plotted in Fig. O Only for q = 1 the 
entanglement entropy is maximal at a stationary point. 



be strictly decreasing functions of k+ . The limiting pro- 
cedure (|6.2[) selects the value of k+ corresponding to the 
maximal entanglement entropy (the circles in Fig. 1201) . 
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